Hopping and microscopic dynamics of ultrasoft particles in cluster crystals 
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We have investigated the slow dynamics of ultrasoft particles in crystalline cluster phases, where 
point particles interact through the generalized exponential potential u{r) = eexp[—{r /a)"], focus- 
ing on the cluster fee phase of this model with n — 4. In an effort to elucidate how the mechanisms 
of mass transport depend on the microscopic dynamics and in order to mimic a realistic scenario 
in a related experiment we have performed molecular dynamics, Brownian dynamics, and Monte 
Carlo simulations. In molecular dynamics simulations the diffusion of particles proceeds through 
long-range jumps, guided by strong correlations in the momentum direction. In Monte Carlo and 
Brownian dynamics simulations jump events are short-ranged, reflecting the purely configurational 
nature of the dynamics. In contrast to what was found in models of glass-forming liquids, the effect 
of Newtonian and stochastic microscopic dynamics on the long-time relaxation cannot be accounted 
for by a temperature-independent rescaling of the time units. From the obvious qualitative dis- 
crepancies in the short time behavior between the three simulation methods and the non-trivial 
difference in jump length distributions, long time relaxation, and dynamic heterogeneity, we learn 
that a more complex modeling of the dynamics in realistic systems of ultrasoft colloids is required. 

PACS numbers: 61.20.Ja, 64.70.pv, 82.70.Dd 



I. INTRODUCTION 

Although the formation of stable clusters of particles 
that interact via entirely repulsive, ultrasoft potentials 
might seem counterintuitive at first sight, it is by now 
well established that cluster phases of such systems do 
exist [iHEl- This phenomenon is potentially relevant for 
a wide class of ultrasoft colloidal systems, such as den- 
drimers and microgels, in which the centers of mass of 
the colloidal particles can overlap with a finite energy 
cost. Theoretical considerations |7| in combination with 
computer simulations Q for models of ultrasoft colloidal 
particles have given evidence for the existence of these 
clusters and have, in addition, provided a deeper under- 
standing of this phenomenon. Stable clusters form either 
in disordered, liquid-like phases or occur as ordered clus- 
ter crystals where the particle aggregates populate the 
lattice sites of regular fee or bee lattices. Most of the 
investigations have been dedicated up to now to study 
the static properties of the cluster phases, of which we 
mention the two most remarkable ones: (i) freezing and 
melting lines depend in a linear way on the density and 
(ii) the lattice constant of cluster crystals is invariant un- 
der compression, inducing thereby a linear growth of the 
cluster size with density. 

In contrast, little effort has been dedicated, so far, to 
obtain a deeper insight into the dynamics that governs 
these cluster forming systems. To the best of our knowl- 
edge, only two contributions that deal with the diffusion 
and the relaxation dynamics in cluster crystals have been 
published d, From these investigations, based on 



molecular dynamics (MD) simulations, it became clear 
that the dynamics in cluster crystals shows features at 
least as intriguing as their static counterparts. For in- 
stance, these investigations revealed that particles move 
constantly between neighboring clusters, while maintain- 
ing the original lattice structure of the cluster crystal and 
the average cluster population on the lattice sites. Fur- 
ther, a closer analysis of the dynamic correlation func- 
tions provided evidence for a decouplingbetween self and 
collective time-dependent correlations [9|. 

The present contribution is dedicated to investigate the 
dynamical properties of ultrasoft particles in a crystalline 
cluster phase in detail and from a more general point of 
view. In particular, we intend to study how mass trans- 
port is realized in different simulation schemes, which 
mimic different scenarios in related realistic systems: (i) 
in MD simulations, where the influence of the microscopic 
solvent is neglected and the mesoscopic (colloidal) parti- 
cles move according to Newton's law; (ii) in Monte Carlo 
(MC) simulations, where the random collisions of the col- 
loids with the solvent particles are mimicked through the 
stochastic nature of the algorithm; (iii) in Brownian dy- 
namics (BD) simulations, where the influence of the sol- 
vent is taken into account via the friction term and the 
stochastic forces acting on the particles. 

The model system that we choose to study cluster 
phases is the generalized exponential model of index n 
(GEM-n) [1], in which particles interact via the poten- 
tial 



(f>(r) = ecxp[-{r/ay 



(1) 
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In the literature, this model has been used to describe 
the effective interactions between linear polymer chains 
{n — 2) [ll| and dendrimers {n « 3) [I4I in the dilute 
regime. For n > 2, we expect this model to capture the 
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general features of purely repulsive, cluster-forming sys- 
tems . In contrast to Ref. Q , where n = 8 was consid- 
ered, we have based our investigations on the GEM-4 sys- 
tem, for which highly accurate data for the phase diagram 
are available 0, Q. In a first step, we have focused on 
the microscopic dynamics and on the jump events, com- 
paring MD and MC dynamics. We have found distinct 
differences between these two types of dynamics imposed 
by the simulation technique: in MC simulations, hop- 
ping takes place only over short distances (i.e., essentially 
from one cluster to a neighboring one); in contrast, in 
MD simulations particles can migrate over large distances 
through the crystal. Furthermore, our investigations on 
the diffusion constant and on the dynamic correlation 
functions have shown that properties calculated from MD 
and MC cannot be superposed by a simple, temperature- 
independent rescaling of time. This finding is in striking 
contrast to what is known from normal [isl [T^ or glass- 
forming 

[iMll liquids. Finally, we have found that BD 
simulations give rise to dynamical properties very simi- 
lar to those observed in MC simulations, supporting the 
view that MC dynamics effectively (although approxi- 
mately) incorporates solvent effects |18l - [20| . Our results 
indicate that a complex dynamical scenario may arise 
as a competition of activated processes and momentum 
correlations, and suggest that the solvent properties may 
change qualitatively the nature of the slow dynamics in 
cluster-forming colloidal systems. 

The paper is organized as follows: in Section |lT] we in- 
troduce the model and present our simulation methods; 
in Section IIIII we present our results for the dynamical 
properties of a GEM-4 cluster crystal and finally in Sec- 
tion IIVI we give our conclusions. In the Appendix we 
present our algorithm for cluster identification. 



II. MODEL AND METHODS 

Our model is composed of classical particles of mass 
m enclosed in a cubic box with periodic boundary con- 
ditions. Particles interact through the generalized expo- 
nential potential, Eq. ([l}, with rt = 4. The potential 
is truncated and shifted at Vc = 2.2a. In the following, 
we will use the parameters a and e in Eq. ([T]) as units 
of distance and energy, respectively. The units of time 
arc chosen differently according to the simulation method 
(see below). For each studied density in the fee cluster 
phase (see Fig. [T]), the equilibrium value of N was deter- 
mined using the algorithm developed in Ref. Q , yielding 
N = 3367 for p = 6.4. During this procedure, we used 
an fee cluster crystal comprising four unit cells per side. 
We have checked that size effects do not play a relevant 
role by performing simulations on smaller systems, i.e., 
using three unit cells per side. 

The dynamical properties of the model were investi- 
gated using MD, MC, and BD simulations. MD simula- 
tions were performed in the NVE ensemble with fixed 
center of mass using the velocity- Vcrlct algorithm [2l| . 



Equilibration at various temperatures was achieved by 
reselecting the velocities of particles at regular time in- 
tervals according to the appropriate Maxwellian distri- 
bution. MC simulations were performed in the NVT 
ensemble using the standard Metropolis algorithm. At- 
tempted moves involved only single particle displace- 
ments, randomly generated over a cube of side 0.6. The 
acceptance ratio ranged from 64 % (at high T) to 42 % 
(at low T). At any fixed T, the acceptance ratio de- 
creases monotonically with the length of the maximum 
attempted displacement. Thus, in contrast to what found 
for a model glass- forming liquid [l5| , it is not possible to 
define an "optimal" value of the maximum attempted 
displacement. The value chosen in this work is a reason- 
able comprise between efficiency of the simulation and 
physical realism. BD simulations were performed by in- 
tegrating a set of coupled, inertia-free Langevin equa- 
tions using a simple Eulcr algorithm [2l|. We used a 
friction parameter ^ = 10^'^ (in reduced units) and a 
temperature-dependent time step (5iBD- ? was chosen 
sufficiently small that the dynamical properties at high 
T no longer depend on ^. To avoid artifacts in the dy- 
namic correlation functions measured during MC and BD 
simulations, the motion of the center of mass of the sys- 
tem was subtracted out from t he part icle displacements. 
In the following, we will use -^/mcr^/e as the unit of time 
for MD simulations. In these units, the time step J^md 
for the integration of the equations of motion was 0.03. 
In the case of MC simulations, time will be measured in 
units of MC sweeps, where one MC sweep corresponds 
to N attempted displacement moves. In BD simulations, 
the time unit will be fi^/Doj where Dq = 1/^ is the diffu- 
sion coefficient at infinite dilution and unit temperature. 
(5tBD was adjusted at each state point so that the dy- 
namical properties did not depend on it anymore upon 
further reduction, and ranged from 5 • 10^® (at low T) to 
10~* (at high T). Procedures to match the time scales 
of the different simulations will be discussed in Sec. IIIII 



III. RESULTS 

In this section we characterize the dynamical proper- 
ties of GEM-4 particles in the fee cluster crystal phase, 
paying particular attention to jump diffusion processes 
and to the effects of the different types of microscopic 
dynamics introduced in Sec. [Til We will focus on the be- 
havior of the system as a function of temperature along 
the isochore p = 6.4 (see Fig. [Ij. If not stated otherwise, 
in the following we will always refer to this particular 
density. The investigated temperature range comprises 
both the stable fee cluster phase (T < 0.80) and the 
super-heated regime. In the latter, the underlying fee 
crystal structure remained (meta)stable throughout the 
simulation. We also performed simulations for selected 
state points at different densities (p = 4.3 and 8.7, see 
Fig.[T|) and checked that the expected scaling of dynamic 
quantities as a function of p/T holds @- In the following 
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FIG. 1: State points studied in this work in the temperature- 
density phase diagram. Triangles indicate state points for 
which the fee cluster crystal was metastable. The lines be- 
tween cluster fee, cluster bcc, and fluid phases of the GEM-4 
system have been redrawn from Ref. 0|. 
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FIG. 3: Time scaling factor fMc/(^MD/'5iMD) as a function 
of pjT . (5tMD is the time step used in MD simulations. t\^Q 
and tJiD are defined as the times at which the mean square 
displacement equals a fixed target value in MC and MD sim- 
ulations, respectively. The target values considered are 4 
(squares) or 25 (circles). 
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FIG. 2: Main panel: mean square displacement {5r'^{t)) as a 
function oft a,t p — 6.4 and T — 0.60 from MD (solid line) and 
MC (dashed line) simulations. Inset: same as main panel but 
as a function of the rescaled time t/t* (see text for definition). 

wc will therefore use p/T as the natural control parame- 
ter of the system. In Sec. IIII Al and IIII Bl we will foeus 
our attention on the comparison between MD and MC 
simulations. In See. IIII CI we will discuss the results of 
BD simulations. 



A. Microscopic dynamics and jump events 

We start our discussion with the analysis of the mean 
square displacement {Sr'^{t)) = (|r(i) — r(O)p). This al- 
lows us to review some general features of the dynam- 
ics discussed previously for GEM-8 particles and to 
perform a first comparison of MD and MC simulations. 
In the main panel of Fig. [2] we show the mean square 



displacement for a representative state point well within 
the stable fee phase (p = 6.4, T = 0.60). After the 
ballistic time range (not shown here), the MD dataset 
shows distinct oscillations due to single-particle vibra- 
tional modes [1, [^, HH, which disappear for t > 10. At 
long times, normal diffusion is recovered, i.e., ((5r^(t)) ~ 
6Dt. Due to the arbitrary choice of the "time unit" in the 
case of MC simulations, the respective curve is shifted. 
It is possible, however, to rescale the time variable, t/t*, 
so that the long time behavior of {5r'^{t/t*)) coincides 
for both types of dynamics. The state-dependent scal- 
ing times t* are defined, for both dynamics, by the con- 
dition (Sr'^{t*)) = 25. The precise choice of the "tar- 
get" mean square displacement is irrelevant as long as 
this value is within the diffusive regime (see also Fig. [3|) . 
The original and the rescaled datasets are shown in the 
main panel and inset of Fig. [21 respectively. Adopting 
the rescaled time representation, the mean square dis- 
placements obtained from MD and MC dynamics nicely 
collapse onto a unique curve at long times {t/t* > 10~^). 
Clear discrepancies are found at short times, due to the 
expected phonon suppression in MC dynamics [l^ . Sim- 
ilar rescaling approaches were used in previous works on 
Lennard- Jones liquids [l3l [T5|. However, in contrast to 
what was found in those studies, the ratio of the MD 
and MC scaling times displays a systematic temperature 
dependence. This is demonstrated in Fig. [31 where we 
show the time scaling factor t^j^^ / (t^^j^ / Stun) for two 
values of the target mean square displacement. The ra- 
tio iMc/(*MD/'^*MD) decreases monotonically with 1/T 
within the explored temperature range. Therefore, it is 
not possible to set, once and for all, the time conversion 
factor between the different types of dynamics. 

To understand the discrepancy between the thermal 
behaviors of the system in MD and MC simulations, we 
analyze in more detail both types of dynamics and iden- 
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FIG. 4: Distribution of net jump length Vnet as a function of 
r/dniL at p = 6.4 and T = 0.60 for MD (empty circles) and 
MC (filled circles) simulations. 



tify the elementary hopping processes leading to particle 
diffusion. To this end, we monitor individual particle 
trajectories and identify jump events by means of the 
cluster analysis outlined in the Appendix. Our definition 
of jump events is as follows. At any instant of time t, 
we map the position r(t) of a tagged particle to the cen- 
ter of mass (CM) of its respective, initial cluster, located 
at RJ^^tiai R""(t). In MD simulations, a particle is 
considered "equilibrated" in its respective cluster if the 
residence time exceeds a characteristic time, t*q, of the 
order of a few vibrational periods (cf. Fig. [5]). In our 
calculations at p — 6.4. we chose t*^ = 3.6, independent 
of temperature. We checked that reasonable variations 
of t*q around the selected value by some 10% or 20%, or 
introduction of a T-dependent t*^ do not modify qualita- 
tively our analysis. In MC simulations, the equilibration 
time is given by ieql^Mc/^MD)' where t^c ^md are 
the scaling times introduced above. A jump event starts 
when the tagged particle leaves the cluster in which it was 
originally equilibrated and ends when the residence time 
of the tagged particle in any of the visited clusters exceeds 
t*q. The net jump vector is then Tnet = R-lnai - R-ii"tiab 
where R^^ai t^*^ position of the final cluster site 
visited during the jump process. This analysis is carried 
out for all jumping particles in the system and averages 
are performed over all identified jump events. 

The normalized distribution of net jump lengths 
Vnctir) = T'drnotI) is shown in Fig. |4]for the same state 
point as considered in Fig. [51 The MD results display 
marked peaks corresponding to integer multiples of the 
nearest neighbor positions in the fee lattice, suggesting 
preferential motion along straight paths. Strikingly, the 
distribution extends up to and even beyond lOdnn, where 
dnn is the nearest neighbor distance. That is, parti- 
cle diffusion in MD simulations proceeds through corre- 
lated, long range jumps, whose typical length exceeds 
(inn- The estimated probability of long-range jumps, 
drVnctif), is approximately 50% at the temperature 



FIG. 5: Main panel: distribution of net jump length Vnct 
as a function of r/dnn for MD simulations at various tem- 
peratures. The dashed line indicates an exponential func- 
tion ~ exp(— r/f), with f ~ 1.43dnn. Inset: log-log plot of 
Vnetir/dnn) for T = 0.70. The dash-dotted line is a power 
law decay l/r^+" with a = 2.2 ± 0.2. 

considered in Fig. UJ while the estimated probability of 
jumping to the nearest neighbor site is ^ 30%. These 
values do not depend significantly on temperature (see 
below). Hence, correlated long-range jumps provide the 
leading mechanism for mass transport during MD simu- 
lations. The shape of the distribution of net jump lengths 
changes considerably when the MC dynamics is consid- 
ered. In this case, 7^not(?') decays rapidly beyond r dnn- 
This can be understood in terms of the stochastic nature 
of the MC dynamics and suggests that long range jumps 
in MD simulations are guided by correlations in the mo- 
mentum direction along the jump path. We will discuss 
this point in more detail below. 

For both types of microscopic dynamics, the net jump 
length distribution Vnetir) displays a mild dependence 
on temperature. In Fig. [S] we show Vnctir) for MD sim- 
ulations at temperatures T = 0.80, 0.70, and 0.55 for 
MD simulations. This temperature range corresponds to 
a variation of diffusivity of approximately three decades 
(see Fig. ini further discussed below). In general, Vnetir) 
decays a bit more rapidly as T decreases. However, the 
probability of long range jumps in MD simulations re- 
mains substantial even at the lowest investigated tem- 
perature (T = 0.55). Hence, the qualitative differences 
between MD and MC simulations discussed above with 
reference to Fig. |4] are relevant for the whole investigated 
T range. At large distances (r > 2dnn) the envelope of 
the maxima of Vnot (f) obtained from MD simulations can 
be described reasonably well by an exponential function, 
aexp {—r/f), over approximately two decades in 7'nct('')- 
f is a characteristic jump length, which depends very 
weakly on T and whose value is « 1.43(inn- 

As it can be seen from the inset of Fig. [5l at even 
larger distances the shape of Vnetif) seems to crossover 
to a power law decay, a behavior reminiscent of 

Levy flights [1^ [13] . The estimated power law exponent 
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FIG. 6: Correlation between angles 6 and net jump lengths r 
along jump paths at T = 0.80 from MD simulations. Empty 
circles: angles 9 between successive jumps versus net jump 
length r. Filled circles: correlation between average values 
of angles at fixed jump length and the jump lengths. The 
solid line indicates the corresponding distribution of net jump 

lengths Pnot(r/dnn) (cf. Fig. El). 




FIG. 7: Possible consecutive jumps to nearest neighbors in 
the fee unit cell. The full black and grey circles indicate fee 
lattice sites. Blue arrow: 60°, r = dnn. Red arrow: 90°, 
r — \/2rfnn. Green arrow: 120°, r = ^/3dnn. Orange arrow: 
180°, r = 2dnn. 



(1 + a w 3.2) is slightly higher than the upper bound 
required for standard Levy flights [1^ [l^l ■ Furthermore, 
inspection of the mean square displacement (see Fig. [2]) 
does not show clear signatures of anomalous diffusion — a 
feature commonly associated with Levy flights [l^ — on 
any time scale accessible by our MD simulations. On the 
basis of the current data, it remains unclear whether a 
Levy flights picture would be appropriate for the descrip- 
tion of diffusion in cluster crystals of ultrasoft particles. 
More extensive investigations are required to completely 
settle this issue. 

To better characterize the nature of the long range 
jumps occurring in MD simulations, we now study the 
relationship between the net jump length r and the an- 
gles 9 enclosed by successive steps between cluster sites 



along the jump path. In Fig. |6]we display the measured 
angles 6 versus the corresponding net travelled distance 
r for a representative state point (T = 0.80). For each 
jump of length r we stored n — 1 angles, where n was the 
number of steps along the path. The horizontal bands at 
60°, 90°, 120°, and 180° reflect the possible angles be- 
tween consecutive steps of length dnn in the fee lattice 
(see Fig. [7|). Angles 6* w occur for consecutive back 
and forth steps, while the data at r « are associated 
to jumps, possibly comprising several intermediate steps, 
whose initial and final clusters coincide. 

Fig. [6] shows that the data tend to accumulate at large 
angles, 9 = 180° and 120°, corresponding to jumps to 
second nearest neighbors of the fee lattice along straight, 
or slightly deflected, directions. The filled circles in Fig.[S] 
indicate the average angles at a fixed average net jump 
length. On average, large angles correspond to large 
jump lengths, revealing that particles tend to proceed 
along approximately straight paths, possibly keeping a 
significant correlation in the direction of their momen- 
tum along the jump path. To illustrate the peculiar 
nature of the hopping processes observed in MD simu- 
lations, we show in Fig.|S]the trajectories of selected par- 
ticles at r = 0.80. They display correlated jumps over 
different lengths, ranging from a few nearest neighbor 
distances (Fig. [8]-a) to over ^ 10 nearest neighbors dis- 
tances (Fig.[5]-c). This latter event clearly shows a strong 
directional correlation of the movement, persisting over 
several intermediate steps along the jump path. 



B. Diffusion, heterogeneity, and relaxation 

Naturally, the question arises: how much are transport 
properties affected by the difference in the microscopic 
dynamics? To discuss this point, we compare in Fig. [9] 
the T-dependence of the diffusion coefficient D obtained 
from MD and MC simulations using the standard Ein- 
stein relation, i.e., {5r'^{t)) = QDt. As it is customary 
for systems displaying slow dynamics, the Arrhenius rep- 
resentation is adopted. In the main panel of Fig. [9] the 
diffusion coefficient for MD, I?md, is multiplied by 0.03 
for better comparison with the MC data. The fact that 
it is not possible to present both sets of data on a single 
curve using a simple rescaling of the time units, demon- 
strated in the figure by the different slopes of the two data 
sets, is consistent with the systematic T-dependence of 
the ratio of scaling times ^mc/^md (cf- Fig- El)- This 
finding is particularly relevant in view of previous nu- 
merical investigations on liquids in equilibrium [isl [l^ 
or in the supercooled regime p"5l - [l7| . In the latter three 
studies, MD and Brownian dynamics [iBl or MD and 
MC dynamics [3, were used. It was shown that the 
long time relaxation of glass-forming liquids is indepen- 
dent of the microscopic dynamics. This means that relax- 
ation times, transport coefficients or dynamic correlation 
functions obtained via the different simulation methods 
can be superposed after a T-independcnt rescaling of the 
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FIG. 8: Representative jump events of tagged particles (indicated by tlie colored small interstitial sphere) during MD simulations 
at r = 0.80. Large particles indicate the centers of mass of the visited cluster sites. The cluster sites visited during the jump 
path are highlighted with large colored spheres. The positions of the tagged particle along the trajectory are marked at 
equidistant steps. 
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FIG. 9: Main panel: Arrhenius plot of the diffusion coefficient 
Dmd from MD (empty symbols) and Dmc from MC (filled 
symbols) simulations. Inset: scaled diffusion coefficients D* 
obtained from MD and MC simulations (see text for defini- 
tion). In both panels, dotted lines represent fits to the Ar- 
rhenius law Do exp{—E/T), where E is an effective activation 
energy. The activation energy E equals 1.57p and 1.44/9 for 
fits to -Dmd(T') and Dmc{T), respectively, and 1.39p for fits 
to D* {T) , independent of the microscopic dynamics. 



time units. For the slow dynamics of ultrasoft particles 
in cluster crystals, however, such a rescaling approach is 
hindered due to the prominent role of momentum cor- 
relations. Hence, our results indicate that independence 
of long time relaxation from the microscopic dynamics 
solely holds for systems where the relevant dynamical 
correlations are configurational in nature. 

MD simulations of a GEM-8 cluster crystal 0] have 
shown that the T-dependence of D closely follows the 
Arrhenius law, D = Doexp{—E/T), where E = E{p) is 
an effective activation energy that scales proportionally 



to p [13 . This observation, combined with our analy- 
sis of the jump length distributions (Figs. |l]and[S|), sug- 
gests a way to rationalize the discrepancy between MC 
and MD results for D. This can be achieved by invok- 
ing a thermally activated diffusion process governed by a 
unique activation energy, i.e., independent of the micro- 
scopic dynamics, but based on different distributions of 
elementary jump lengths. Towards this aim, let us write 
down the Arrhenius law for activated diffusion assuming 
a distribution Vnctir) of jump lengths 



D 



exp{~E/T), 



(2) 



where (t^) is an average waiting time, which ultimately 
sets the natural time scale for the process. As long as 
the T-dependence of (tw) is approximately similar for 
the different dynamics, it should be possible to rectify 
the Arrhenius plots of Fig. [9] by considering the rescaled 
diffusion coefficient D* ~ D / J^r'^'Pnctir)dr. The in- 
set of Fig. |9] shows that this is indeed the case. Fits 
to an Arrhenius law for D*(T) provide a unique effec- 
tive activation energy E w 1.39p, independent of micro- 
scopic dynamics. We found that the fitted value of E is 
close to the barrier height Ei, separating neighboring clus- 
ters, estimated [l^] from single particle potential energies 
[Ef, « 0.94i5). This clearly points to a non-cooperative 
nature of activated dynamics in the system under con- 
sideration. 

The microscopic dynamics exerts even a stronger in- 
fluence on other typical time-dependent correlation func- 
tions of the liquid state, such as density-density time cor- 
relators. In Fig. [10] we show the self part of the van 
Hove correlation function Gs{r,t) at T = 0.80 and com- 
pare MD and MC datasets for t = t* , where again the 
scaling time t* is chosen such that (Sr^{t*)) = 25. The 
MC dataset can be simply described by a Gaussian func- 
tion with variance 2Dt (as expected for standard random 
walks) modulated by peaks centered at the nearest neigh- 
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FIG. 10: Main panel: Self part of the van Hove correlation 
function Gs{r,t = t*) as a function of r/dnn at T = 0.80 
for MD (thick solid line) and MC (thick dashed line) simula- 
tions. The time t* is chosen such that {5r'^{t*)) = 25. The 
thin dashed line indicates the expected functional form for 
standard random walk, inr exp[— r ^/{ADt)]/{4TvDtf^^. The 
thin solid line indicates an exponential function ~ exp(— r/f), 
with f ~ 2.44rfnn. Inset: same as main plot but on a semi- 
logarithmic scale. 




FIG. 11: Non-Gaussian parameter 02 as defined in Eq. (|3]) 
as a function of the rescaled time t/t* at T = 0.60 from MD 
(solid line) and MC (dashed line) simulations with t* as de- 
fined in the text. 



bor distances of the fee crystal. The MD dataset, on the 
other hand, is much more stretched and displays an ap- 
proximately exponential envelope. This behavior is a di- 
rect consequence of the broad distribution of elementary 
jump lengths (cf. Fig. |4] and [5]) . 

The stretched shape of the van Hove correlation func- 
tion obtained from MD simulations motivates a closer 
inspection of the non-Gaussian parameter 
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FIG. 12: Maximum value of the non-Gaussian parameter 
af'' as a function of p/T for MD (empty circles) and MC 
(filled circles) simulations. 



context of glass- forming liquids, a2(t) has also been 
used l25l . l26ll as a simple measure of "dynamic heterogene- 
ity" [ll-ia- In Fig. Owe show a2{t/t*) at T = 0.60. 
At short times qualitative discrepancies appear between 
MD and MC data, as expected: marked oscillations in 
a2{t/t*), due to single-particle vibrations, are observed 
in fact only for MD simulations. At longer times, on the 
other hand, the overall shape of a2{t/t*) is very similar 
to the one found in normal and supercooled liquids [2^ : 
for both types of simulations, a2 displays a peak at inter- 
mediate times corresponding to a maximum value a^^^. 
As in standard glass-formers, a™'*'^ grows by decreasing 
T, albeit in a moderate fashion (see Fig. [T^ . However, 
the values of al^^^ obtained in MD simulations are signif- 
icantly larger (up to a factor 40) than in the MC simula- 
tions of our system or than the typical values observed in 
MD simulations of model glass- forming liquids [2^ . This 
is a consequence of the broad distribution of elementary 
jump lengths found in MD simulations for the system 
under consideration. Such a dynamic heterogeneity is 
essentially non-cooperative in origin and therefore differ- 
ent from the one normally observed in glass-forming liq- 
uids [27l - l29j . By visual inspection of animated trajecto- 
ries, we have checked that jump events are not associated 
to correlated motions of several particles (as observed, 
for instance, in [l^l), but rather reflect the migration of 
individual particles over the sample. 

Finally, we study the dependence of the self interme- 
diate scattering function 



1 ^ 

= ]^ E ( '^^P ■ [r, it) - r, (0)] } ) 



which quantifies the deviation from Gaussianity of the 
distribution of particle displacements at time t. In the 



on the microscopic dynamics. Due to the underlying 
crystal structure, the relaxation of Fs(k,t) is strongly 
anisotropic, i.e., it depends explicitly on the direction of 
k. While this effect of anisotropy is expected to be most 
pronounced at reciprocal lattice wave-vectors, we found 
in addition that Fs{\i,t) possesses a non-trivial angular 
dependence at any fixed norm (not shown here) . To sim- 
plify the discussion, in the following we focus only on 
directions in reciprocal space corresponding to the first 
and second reciprocal lattice wave-vectors of the fee lat- 
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MC 



MD 




FIG. 13: Self intermediate scattering function Fs(k,t) as a function of t at T = 0.80. The wave- vectors considered are chosen 
along the directions ki (panels (a) and (b)) and k2 (panels (c) and (d)) corresponding to the first and second reciprocal lattice 
wave- vectors of the fee lattice, respectively. The norm of k is varied from to 2ki, with i = 1,2 (curves from bottom to top); 
the curves are colored according to the colorscale shown on the left side of the figure. Left and right panels are results from 
MC and MD simulations, respectively. 



tice, ki/Zci and k2/fc2. Along each of these directions, we 
vary the norm of k within the range < k < 2ki , where 
ki is the norm of the selected reciprocal lattice vector 
(fci = 5.39, k2 = 6.23). Results obtained at T = 0.80 are 
collected in Fig. [131 Clearly, the relaxation of (k, t) be- 
comes extremely slow at reciprocal lattice vectors. How- 
ever, the relaxation patterns for off-lattice wave-vectors 
are rather complex in the two dynamics: as the norm 
of k approaches that of a reciprocal wave-vector the 
relaxation slows down, but the difference in relaxation 
times between "fast" and "slow" wave-vectors depends 
markedly on the microscopic dynamics. Our results in- 
dicate a larger heterogeneity of relaxation times for MC 
dynamics, a fact which contrasts our previous analysis 
of the non-Gaussian parameter. A similar conclusion is 
drawn from the analysis of the angular dependence of 
-F's(k, t) at fixed norm of k (not shown here). Further 
investigations are thus required to clarify this behavior 
and its possible connection to the decoupling between self 
and collective relaxation discussed in Ref. M . 



C. Comparison with Brownian dynamics 

Our analysis has focused so far on the comparison of 
dynamical properties obtained through MD and MC sim- 
ulations. MC simulations are expected to mimic in some 
sense the Brownian motion of colloidal particles due to 
collisions with solvent particles. Simulation data [l^ and 
approximate analytical models [20j provide support for 
this view, at least for small MC displacements. In this 
subsection we explicitly compare BD and MC simulations 
for the GEM-4 model, and show that the two methods 
lead indeed to very similar macroscopic dynamics. 

We start with a direct assessment of the net jump 
length distribution Vnct{i" / dnn) obtained from BD simu- 
lations at r = 0.60. The results are displayed in the main 
panel of Fig. I14| and compared with the corresponding 
dataset from MC simulations already presented in Fig. 2] 
The net jump length distribution from BD simulations is 
short-ranged and remarkably similar to that obtained us- 
ing MC simulations. A small peak at r « 2dnn indicates 
rare jumps over two nearest neighbor distances, a feature 
absent in the MC dataset. In the inset of Fig.[T4]we show 
the van Hove correlation function Gs{r, t*) calculated at 
a time t* at which the MSD equals 9. Wc observe a 
striking correspondence between BD and MC datasets, 
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FIG. 14; Main panel: Distribution of net jump length "Pnct as 
a function of r/dnn at p = 6.4 and T = 0.60 for BD (empty 
squares) and MC (filled circles) simulations. Inset: Self part 
of the van Hove correlation function Ga(7", t — t*) as a function 
of r/dnn at T = 0.60 for BD (solid line) and MC (dash-dotted 
line) simulations. The time t* is chosen such that {Sr^{t*)) = 
9. 
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FIG. 15: Main panel: Arrhenius plot of the diffusion coeffi- 
cient 0.003 • Dbu from BD (empty squares) and Dmc from 
MC (filled circles) simulations. The scaling factor for Dbd 
has been chosen so as to maximize the overlap of the two 
data sets. Inset: time scaling factors fMc/(^BD/0.003) (open 
squares) and t^^c / {'^Iau / ^'^md) (filled circles) as a function of 
p/T. fJiCi ^mdj BLiid t^Y) are defined as the times at which 
the mean square displacement equals 4 in MC, MD, and BD 
simulations, respectively. 



demonstrating a close similarity between the dynamics 
emerging from the two simulation methods. 

The temperature dependence of the diffusion coeffi- 
cient D{T) is analyzed in Fig. [T5] employing the Ar- 
rhenius representation. Within the temperature regime 
explored by BD simulations, the T-dependences of the 
diffusion coefficients obtained from BD and MC simula- 
tions closely follow one another, up to a trivial, state- 
independent rescaling of the time unit. Thus, not only 
the qualitative aspects of transport mechanisms, but also 



the T-dependences of transport coefficients generated by 
BD and MC are very similar. To contrast the results 
obtained through these stochastic methods to those of 
Newtonian dynamics, we report in the inset of Fig. fTSl the 
ratio of the scaling times ^mc/^bd ^'^d ^mc/^mdj de- 
fined in Sec. IIII M As expected, iMc/^BD remains fairly 
constant throughout the available temperature regime. 
Within the same T-range, ^mc/^'MD displays a system- 
atic variation, as already evidenced in the main panel of 
Fig.im 

Overall, the results shown above indicate a strong sim- 
ilarity between the transport mechanisms at play in BD 
and MC simulations and corroborate the analysis out- 
lined in the previous subsections. In this respect, the re- 
sults shown in Fig.[T4land[T5lsupport the use of MC simu- 
lations as an effective way to incorporate solvent effects in 
a computer simulation of model colloids [l^, H^l • We re- 
mark that an attempt to rcscale the time unit in our MC 
simulations by the acceptance ratio, as suggested in [20| . 
resulted in slightly poorer agreement between the two 
datasets of D{T). Thus, more extensive investigations on 
the connection between dynamical properties generated 
by MC and BD simulations and possible system-specific 
aspects are required. 



IV. CONCLUSIONS AND OUTLOOK 

We have performed detailed MD, MC, and BD simu- 
lations for ultrasoft particles that form stable clusters 
of overlapping particles, which, in turn, populate the 
lattice sites of a regular fee lattice. Our investigations 
on the dynamical properties of the system (in terms of 
particle moves, diffusion and dynamic correlation func- 
tions) revealed striking differences between Newtonian 
and stochastic simulation methods. 

(i) In MD simulations, particle diffusion is realized as 
long-ranged sequences of jumps from one cluster to a 
neighboring one, assisted by strong correlations in mo- 
mentum directions. The envelope of the jump length 
distribution is within a good approximation an exponen- 
tially decaying function of the distance and extends up to 
and beyond ten nearest neighbor distances. At large dis- 
tances, the jump length distribution is described rather 
well by a power law, l/r^+", with a « 2.2. In MC 
simulations, by contrast, hopping events are limited to 
the nearest clusters, as a result of the purely configura- 
tional nature of the dynamics, (ii) The diffusion con- 
stant, D, evaluated in MD and MC simulations cannot 
be matched on a single master curve by a simple rescal- 
ing of the respective time units. This finding is in strik- 
ing contrast to related results obtained in investigations 
on equilibrated [3, and supercooled [TB - flTj liquids. 
However, by normalizing D with respect to the under- 
lying jump length distribution, an activation energy can 
be filtered out that is common to both types of dynam- 
ics, (iii) Also for the dynamic correlation functions we 
find distinct differences, which become apparent in the 
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overall shape of the self part of the van Hove correlation 
function, in the non-Gaussian parameter, and in the long 
time decay of the self intermediate scattering function for 
off-lattice wave vectors, (iv) BD simulations give rise to 
dynamical properties very similar to those observed dur- 
ing MC simulations, supporting the view that solvent 
effects arc implicitly and effectively included in MC sim- 
ulations through the stochastic nature of the algorithm. 

The present contribution — along with Refs. [1, [l^] — 
can only be considered as a first step towards a deeper 
understanding of the evidently complex dynamics of clus- 
ter forming, ultrasoft systems. Quite a few issues re- 
main unaddressed, which will be dealt with in future 
contributions. First, a more systematic study of the 
influence of the solvent, using full Langevin dynamics 
and possibly including hydrodynamic interactions, is re- 
quired and will be performed in future work. Other 
simulation methods, such as dissipative particle dynam- 
ics [30I Isij . might provide a different, complementary 
view on our data. Another aspect is related to the un- 
derlying model: our inter-particle potential represents an 
effective interaction between two mesoscopic (colloidal) 
particles, where the degrees of freedom of the microscopic 
constituents of the colloids have been traced out (33 ]. 
An important issue is thus related to the level of coarse- 
graining, ranging from the present picture of an "effec- 
tive" particle [12, , over partially averaging over sub- 
entities, such as in multi-blob representations [1^, or 
force- matching schemes [13, H^, to the fully atomistic 
level where all the constituent entities of the colloidal 
particles and the particles of the solvent are considered 
explicitly. Another fundamental question addressing the 
interpretation of the dynamic correlation functions is how 
the dynamics of a system of "effective" particles is related 
to that of a system of atomistically resolved particles. 

In conclusion, our work has revealed qualitative differ- 
ences in the dynamics of ultrasoft, cluster-forming par- 
ticles obtained using Newtonian and two stochastic (i.e., 
Brownian dynamics and Monte Carlo) simulation meth- 
ods. This effect has been attributed to the competition 
between activated processes and solvent effects, which 
controls the distribution of jump lengths and ultimately 
the dynamic correlation functions. We speculate that 
the interplay between activated slow dynamics and sol- 
vent effects could generate complex dynamical scenarios 
in ultrasoft colloids, such as dendrimers and microgels, 
and even a broader class of cluster-forming colloids, in- 
cluding systems with competing interactions. With this 
contribution we hope to motivate further numerical and 
experimental studies on the dynamics of such colloidal 
systems to test this hypothesis. 



p = 4.3. T = 0.75, liq 
p = 4.3. T = 0.50, fee 
/9 = 6.4. T = 0.80, fee 
p = 6.4.T= 1.05, fee 




FIG. 16: Radial distribution function, g{r), of a GEM-4 sys- 
tem as a function of r for different densities p and tempera- 
tures r, both in the liquid and fee cluster crystal phase. From 
the inset we see that g{r) does not vanish completely at its 
first minimum, rmin ~ 0.75. 
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Appendix A: Cluster identification algorithm 

In an effort to trace the migration of the particles 
through the cluster crystal we developed an algorithm to 
distinguish between different clusters and to identify the 
affiliation of a particular particle to a cluster in an unam- 
biguous way at every step of the simulation. In previous 
investigations on clustering systems [3], the position of 
the first minimum in the pair distribution function g{r), 
Tinim was used to provide information whether two parti- 
cles, separated by a distance r, belong to the same cluster 
(for r < Tniin) or not (for r > rnii„). This criterion can 
undoubtedly be used to obtain a first, rough estimate 
for identifying those particles that belong to a particu- 
lar cluster. However, a closer analysis of g{r) reveals (cf. 
Fig. I16p that even for relatively low temperatures this 
function does not vanish completely for r ~ r,^in- Par- 
ticles migrating between two neighboring cluster sites of 
the cluster crystals are to be made responsible for this 
effect. 

In the following we present a refined version of the 
cluster identification algorithm presented in Ref. 0; it 
can be subdivided into the following four steps: 
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We start our procedure with the first particle and 
identify all particles separated by a distance less 
than a given cut-off radius rc as neighboring parti- 
cles. This procedure is repeated for all remaining 
particles. As a first guess we choose rc = rmin, a 
value which will be corrected iteratively during this 
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FIG. 17: Procedure to separate merged clusters, (a) Three 
neighboring clusters are merged by hopping particles. At this 
step, all particles in red form a single cluster, (b) The search 
for particle neighbors is repeated for the merged cluster with 
a reduced cut-off radius Vc. (c) Particles (light grey spheres) 
with a small number of neighbors {ric < n™'") are excluded, 
(d) The search for particle neighbors is repeated ignoring the 
light grey spheres, identifying thereby disjointed clusters, (e) 
Excluded particles (light grey spheres) are reintroduced and 
reassigned to the respective nearest clusters, (f) Procedure 
accomplished, yielding three disjoint clusters. 



algorithm. 

With this information at hand, we start again with 
the first particle and label all its neighbors, their 
respective neighbors, etc.; in this way, we are able 
to identify a first cluster of the system. We then 
proceed to the next particle that has not been la- 
belled yet, and repeat the operation. Finally, all 
particles have been assigned to a cluster. 

At this stage, the algorithm reproduces exactly the 
results obtained in [3]. It risks, however, to provide 
misleading data: as particles move from one cluster 
to another, the particles of these two clusters might 
now be counted to belong to the set of neighbors 
of the hopping particle, merging thereby the two 
clusters (see Fig. [TTl-a). In our refined algorithm, 
we modify the first guess of Tc in an iterative way 
and eventually all merged clusters and all single- 
particle clusters are eliminated. To this end, we 
reduce Tc and introduce three check parameters: 
n™'" and n^^^, the expected minimum and max- 
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FIG. 18: Distribution of cluster sizes for various state 
points. Upper panel: distribution calculated after Step 2 
of the cluster analysis algorithm. The inset shows an en- 
largement of the large-ric region, revealing the occurrence of 
merged clusters. The arrow indicates the occurrence of clus- 
ters composed by few particles. Lower panel: distribution cal- 
culated after Step 4. The inset reveals the absence of merged 
clusters. The arrow indicates that all clusters of size less than 
n!?"" have vanished. 



imum cluster size present in our system (roughly 
estimated from the cluster size distribution calcu- 
lated after Steps 1 and 2) and Nc, the number of 
lattice sites in our system which is, by definition, 
equal to the amount of clusters in the crystal phase. 

We proceed as follows: all clusters that have just 
been identified are re-considered. If the size of one 
of them exceeds n^^^, the corresponding collection 
of particles is isolated from the rest of the sys- 
tem (see Fig. [iTl-a) and the search for neighbors 
is repeated within the remaining set of particles 
(Fig. 1171 -b). In a similar manner, those particles 
with the lowest number of neighbors are also iso- 
lated from the others (Fig. [TTI-c) since these parti- 
cles are found to be responsible for merging neigh- 
boring clusters. Ignoring the excluded particles for 
a moment, the neighbors are identified once again 
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as described above (cf. Step 2), giving disjoint clus- 
ters (Fig HTl-d) . Now the isolated particles are re- 
integrated into the ensemble and are reassigned to 
those clusters with the nearest center of mass po- 
sition (Fig, \l7\e and flTl-f). Finally, the following 
checks are made: (i) Does the size of all newly es- 
tablished clusters lie within n™^^ and n™™? (ii) 
Is the number of clusters identified in the system 
equal to the number of lattice sites, Nc? If one of 
these two conditions is violated, the procedure is 
iterated reducing the cut-off radius at each iter- 
ation step. 

4. In a final check on the size of the clusters, those 
particles or small collections of particles that pos- 
sibly have not been assigned to none of the clusters 
before, are assigned to the cluster with the nearest 
center of mass. The success of the improved clus- 
ter analysis algorithm becomes obvious from the 



ensuing cluster population distribution. Peaks due 
to single particles and merged clusters have van- 
ished (see Fig. [HI right panel) , reflecting the cor- 
rect analysis of the cluster sizes distribution of the 
system. 



Finally, to match the identity of clusters at different 
times, we exploit the fact that in the crystal phase the 
centers of mass of the clusters are fixed to their respective 
lattice sites. To be more specific, we found that the root 
mean square displacement of the centers of mass does 
not exceed ten percent of the nearest-neighbor distance 
between clusters, dnn, at any state point in the fee region 
of the phase diagram. We can thus keep track of the 
identity of any selected cluster by locating the cluster 
that, at the next time step, has the closest center of mass 
and by repeating this operation for all time steps. 
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